An extension of RPA preserving energy weighted sum rules. 
An application to a 3-level Lipkin model. 



M. Grasso and F. Catara 
Dipartimento di Fisica, Universita di Catania 



anc 



Istituto Nazionale di Fisica Nucleare, Sez. di Catania 
Cor so Italia 57, 1-95129 Catania, Italy 

Abstract 

A limitation common to all extensions of RPA including only particle-hole 
configurations is that they violate to some extent the Energy Weighted Sum 
Rules. Considering one such extension, the improved RPA (I RPA), already 
used to study the electronic properties of metallic clusters, we show how it 
can be generalized in order to eliminate this drawback. This is achieved 
by enlarging the configuration space, including also elementary excitations 
corresponding to the annihilation of a particle (hole) and the creation of 
another particle (hole) on the correlated ground state. The approach is tested 
within a solvable 3-level model. 



I. INTRODUCTION 

Collective excitations are a common feature of a large variety of many body systems. 
Their properties are intimately related to the stucture of the ground state upon which they 
are built. The simplest theory of excited states of a quantum system where correlations 
are taken into account to some extent is the Random Phase Approximation (RPA). In 
this theory one introduces a set of operators Ql, whose action on the ground state l^o) 
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creates the collective excitations, while the ground state itself is the vacuum for the Q u 
operators. The latters are defined as linear superpositions of particle-hole (ph) creation and 
annihilation operators, occupied (h) and unoccupied (p) single particle states being defined 
with respect to the Hartree-Fock ground state \HF). The X and Y coefficients of these 
linear forms are solutions of equations which can be derived by using the equations of motion 
method [|],|]|. If the hamiltonian contains one- and two-body terms, the solution of these 
equations would imply the evaluation of one- and two-body density matrices. Standard RPA 
is obtained by replacing them by those calculated in the uncorrelated ground state \HF). 
This approximation introduces a visible inconsistency since, on one hand, the definition of 
|^o) as the vacuum of the Q u operators is used to derive the formal equations determining the 
X and Y amplitudes; while, on the other hand, \HF) is used instead of \^q) i n calculating 
the expectation values appearing in those equations. 

Various attempts have been made to eliminate this inconsistency. We quote the pio- 
neering works H|J where the renormalized RPA (RRPA) was introduced. The RRPA was 
applied to study the low-lying spectrum and the transition densities of vibrational nuclei 
H and the double beta decay more recently. A very important contribution to the 
solution of this problem has been given in ||, where a general scheme, the self consistent 
RPA (SCRPA), was developed (see also [|l(J and references therein). In reference f| a 



fully renormalized RPA (fully RRPA) has been proposed, which shares some similarities 



with the approach we are going to present in this paper. In [11,12] it was shown that by 



using the number operator method [13| it is possible to get for the X and Y coefficients 
a closed set of equations having the same form as in RPA, where the density matrices in 
the correlated ground state are expressed in terms of the X and Y coefficients themselves. 
Thus the equations to solve are non linear and their solution requires a big computational 
effort. In order to appreciate how much a better treatment of correlations modifies the RPA 
results, in the same paper a simplified version of the approach was proposed, the IRPA, 
based on the linearization of the equations of motion. The simplification consists in con- 
tracting the two-body terms appearing in the commutator of the hamiltonian with a one 



body ph operator with respect to the correlated ground state. In this way only one body 
density matrices have to be evaluated: the so obtained equations are still non linear, but 
they are much easier to be solved, since only the one-body density matrix appears. When 
the latter is calculated in \HF) rather than in the correlated l^o); RPA is again obtained 
This approach was applied in ]rT| , |T2"|] to the study of the electronic properties of some 



simple metal clusters, obtaining a better description than RPA. However, the formulation 
is quite general and its applicability is by no means limited to such systems. 

A limitation common to all extensions of RPA including only ph configurations is that 
they violate to some extent the Energy Weighted Sum Rules (EWSR). In the present 
paper we will show that this drawback can be eliminated by enlarging the configuration 
space, including also those configurations corresponding to the annihilation of a particle 
(hole) and the creation of another particle (hole) on the ground state. This is in the same 



spirit of Hi5yT6| , where, for the first time, the particle-particle and hole- hole configurations 
were included within the SCRPA approximation. 

Very recently a paper |T7[ came to our knowledge, where the same problem is tackled 



and studied within a solvable 4-level model with a separable residual interaction. As we will 
show below, there are several differences with the present paper: 

1) we explicitly show that the EWSR is exactly satisfied when the configuration space 
is enlarged; 

2) by comparison with the exact solutions of the model we can judge about the quality 
of the results obtained in IRPA and its enlarged version, with respect to the RPA ones; 

3) this comparison allows us to point out that, besides the merit of solving the EWSR 
problem, the approach has the shortcoming that spurious solutions appear. This problem 
is not discussed in |T7] where, indeed, probably because a separable residual interaction is 



used, only one collective state is found despite the fact that 3 elementary excitation modes 
are present in the model. In this context, it is worth mentioning that spurious solutions are 
also found in ||, where they are interpreted as "new excitation modes". 

The paper is organized as follows. In Section II we shortly recall the main IRPA equa- 



tions, pointing to the origin of the EWSR violations. Then we show how this problem is 
solved when the enlarged space is considered. In Section III we illustrate the approach by 



applying it to a solvable 3- level model [IS -21] and comparing the different approximations 



among themselves and with the exact results. 



II. FORMULATION OF THE APPROACH 



In this section we recall the main steps leading to the IRPA equations, presented in detail 



m 



liyi2|| , and illustrate why, being limited to ph excitations, the IRPA approximation 
violates the EWSR 0. Then we show that, enlarging the space by including also pp and 
hh configurations, this difficulty is overcome. In this respect, our approach is similar to the 
fully RRPA [§. 

A. IRPA and the EWSR problem 

Let \^ ) be the ground state of the system and \^ u ) its excited states. Assuming that 
the latters are linear combinations of ph and hp configurations built upon \ty ) one writes: 

|¥„) = Qt|* ) = Y\ X ph B U - W/JI*o> > (1) 

ph 

where p (h) denotes the quantum numbers of an unoccupied (p) and occupied (h) single 
particle state in the uncorrelated Hartree-Fock reference state \HF). In eq.(|l|) we have 
introduced renormalized ph creation (B') and annihilation (B) operators. In ]Tl]JT2] it is 
shown that in the basis diagonalizing the one-body density matrix they can be written as: 

B\ h = D-Ja\a h , (2) 

with: 

D ph = n h -n p , (3) 

where nh and n p are respectively the hole and particle occupation numbers in the correlated 
ground state l^o)- Assuming that l^o) is the vacuum of the Q u operators, 



Qv\*o) = 



(4) 



the ortonormality conditions for the excited states leads to: 



(5) 



The equations determining the X u and Y u amplitudes and the excitation energies E v of the 
states \fy v ) are obtained by using the equations of motion method Jl]||. They read: 



( A B\ 



\ 



B* A* 

with the A and B matrices given by: 



\ Y "/ 



E„ 



(6) 



Aph, P 'h' = (^o\[B p h, H, B p , h ,\^ ) 



and 



B 



ph,p'h' 



-(%\[Bl h ,H,Bl, hl ]\* ) . 



(7) 



(8) 



In (0) and (|8|) H is the hamiltonian of the system and: 



[A,B,C] = -([A, [B,C]} + [[A,B},C]) 



(9) 



The standard RPA equations can be obtained by putting n,h=l, n p =0 in the expressions for 
the operators B and (|2|) and by replacing the correlated ground state appearing in 
(^) and (H) with the Hartree-Fock one \HF). In [[L4[] it is shown that the RPA equations 
can equivalently be obtained by linearizing the commutator [H, B p , h ,] in (|7|) and (H), i-6. by 
contracting it with respect to \HF). A better approximation is done in I RPA, where the 
linearization is made by contraction in I n a loose notation, this means: 



[H, alah] — > a^a + a) a) 



aa 



a f a + ($ \a}a\y }a)a . 



(10) 



Therefore, the occupation numbers in the correlated ground state appear in the IRPA 
expressions, while those in \HF) (i.e. or 1) appear in standard RPA. This procedure 
leads to: 



i _ 1 z' n 1 / 2 n- 1 / 2 i n 1 / 2 n- 1 / 2 ^ \ \ 

Aph,p'h' — T^K^ph U v 'h' + U p'h' U ph ){ e p'pOhh> - thh'O; 



vv 



+D 1 p / h 2 D 1 JZ(hp'\H 2 \ph') , (11) 



where: 



e p , p = {pWp) + Y,n a (p'a\H 2 \pa) (12) 



and 



e w = (hlH^h') + Y / n a (ah\H 2 \ah') . (13) 

a 

For the matrix B one gets: 

B*j W = D^D^hh'^pp') . (14) 

In the above equations H\ is the one-body term of the hamiltonian and H 2 its two-body 
part. We denote by a a generic single particle state (occupied or unoccupied in \HF)). 
As shown in Appendix A of [IT], using the number operator method jIJ3] the occupation 



numbers appearing in the A and B matrices can be expressed in terms of the X and Y 
amplitudes as: 

n P = XK^' ~ 2 ^ D Pi h i X Pih 1 X p*h 1 ) D ph Y p U h Y ph* > (I 5 ) 
/if!/' Pl/ll 

n h = l- £(<W - ^ E ^^A^X^J^F;,^'* . (16) 

pvv 1 Pihi 

Therefore eq.s @ are non linear. They have been solved iteratively in the case of metallic 



clusters |1T| , |I2| . It is, however, apparent that the approach is quite general and can be 
applied to any many body system. The matrices A and B in IRPA, eq.s(|ITj) and ([14]), are 
different from those in standard RPA. On one side the Hartree-Fock single particle energies 
appearing in the A matrix of RPA are replaced by the quantities appearing in the first line 



of eq. QTT|) . On the other side, the residual interaction in the expressions for A and B is now 
renormalized by the factors D 1 ' 2 ^. In RRPA only the latter modification is present. This 
latter modification is present also in RRPA. 



A serious problem arises with respect to the EWSR. As it is well known, if and 
\^ u ) are a complete set of exact eigenstates of the hamiltonian, with eigenvalues E and E u , 
the following identity holds: 

J2(Eu-E )\(%\F\^ )\ 2 = ^ \[F,[H,F]]\^ ) , (17) 

where F is any hermitian single particle operator. The equality (|T~Tf) is in general violated to 
some extent when \^o), and E v are calculated within some approximation. To which 
extent it is satisfied is a measure of the adequacy of the approximation. A very important 
feature of RPA is that eq. (|17]) is satisfied for any one-body operator if, in calculating its two 
sides, one considers \HF) instead of |^ ) an d the solutions of RPA for and [E v — E ) 
|T8| . This feature follows from the fact that, when \HF) is used in (|17D instead of |^ ) 
only particle-hole matrix elements remain in the r.h.s. It is easy to show [[□]] that, if the 
transition operator F has only p — h matrix elements, the two sides of eq.(|17D are equal also 
within IRPA. However, this is not the case in general. 

Let us consider separately the two sides of eq.(|I7|) in IRPA, with a general one-body 
hermitian operator F: 

F = Y,U a « a f>- ( 18 ) 

a/3 

The l.h.s. is easily calculated and gives: 

J2(K - E )\(t> u \F\t> )\ 2 = - £ )K*o|Q,W| 2 

V V 

= J2(E„-E )\(y \[Q u ,F]\y )\ 2 

V 

= - e )\ Y.UhD l £{x; h + r; h )\ 2 , (19) 

v ph 

which is formally equal to the RPA result, apart from the factor D h . Therefore, only the 
ph components of F enter. This is due to the fact that the excited states are described as 
superpositions of ph configurations only. Starting from eq.s(^) and using the properties of 
the X and Y amplitudes, eq.(|T9"D can be written as: 
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^2(E u -E )\(^ u \F\^o)\ 2 

V 

— E D lhfph E D l J'h'fp'h'{A p h, P >h> - B ph y h >) . (20) 

ph p'h' 

In order to evaluate the r.h.s. of eq.flT7|) one can use for the commutator [H, F] the same 
linearization procedure already used in deriving eq.s(|ll|) and (|T4]). It is easy to realize that 
the result of such calculation cannot be equal to ( p0|) since not only the ph matrix elements 
of the residual interaction will appear in it, but also other terms if they are present in the one 
body operator F. This happens because the expectation value of the double commutator is 
taken in the correlated ground state \^o)- We will show this in the next subsection, where an 
enlarged configuration space, including also pp and hh components, will be used to express 
the excited states. Of course, if the correlations present in |\l/o) are small and the occupation 
numbers do not differ too much from and 1, the violations of the EWSR are small. But, 
in general, this is not the case. For example, for Na clusters, the discrepancy was found 



11.12] to be about 25%. 



B. The enlarged space 

As shown in the previous subsection, the problem of violations of the EWSR arises 
because also in IRPA, as in RPA, the excited states are expressed as superpositions of ph 
configurations. Let us then consider the more general expansion: 

l**> = Ql\*o) = E(Kp B U - Y&Baf,)^) , (21) 

a>/3 

where a and (3 stand for any single particle state and a > [3 means that we order these states 
according to decreasing occupation numbers, i.e. n a < np. The operators and B a p are 
an obvious generalization of eq.s (0) and (^|). As before we define |\I/o) as the vacuum of the 
Q u operators: 

Q,|^o) = 0. (22) 

In order to make simpler the notation, we will omit the bars in the collective operators and 
in the states, which, of course, are different from those considered in IRPA since now pp and 
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hh configurations are included, in addition to the ph ones. The derivation of the equations 
of motion can be done by following the same linearization procedure as before. They have 
the same form as in eq.(|6]), the matrices A and B being now: 

1 / n l/2 n -l/2 , n l/2 n -l/2^ 



) 1/2 D V 

'a/3 



+D 1 J 2 Dl /2 ((3 1 \H 2 \a5) (23) 



and 



+D 1 J p 2 Df(5P\H2\ia) , (24) 



where: 



e Q/3 = (alH^P) +5» 7 |#2|/37) ■ (25) 

7 

Apart from the fact that in eq.sfl23|) and (|24|) the indices run over all single particle states, 
the main difference with eq.s(|ll|) and ([L4]) is the presence of the e terms also in the B matrix. 
Coming back to the EWSR problem, eq.([2fj) is easily generalized to: 

^2(E U -E )\(^ U \F\^ )\ 2 

v 

= f<*P D ai3 Y fisDlj 2 {A a p iJ& - B af3a5 ) , (26) 

a>,8 7><5 

which, after some tedious manipulations, can be written as: 
J2(K-E )\(* V \F\* )\ 2 

V 

= 2 Y f^D af3 Y f-ys[^f3S - + ({3j\H 2 \a5)D lS } . (27) 

a/3 75 

The double commutator is easily calculated by using the same linearization procedure 
adopted to derive the equations of motion. Doing that one realizes that eq.flTTD is in- 
deed satisfied. Thus one obtains a kind of generalization of the Thouless theorem. Namely, 
eq.flTTD is satisfied if one calculates its two sides by using the solutions of the equations of 
motion and by making the same approximations introduced in the derivation of the latters. 
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In principle the new approach does not appear to be more difficult than IRPA and its 
equations can be solved in realistic cases by the same iterative procedure used there. In 
practice, however, the computational effort is much heavier since the configuration space is 



much larger. For this reason we have decided to apply it to a solvable 3-level model [19-21 



We show this application in the next section, where we compare the results of IRPA 
and of its enlarged version with the exact solutions of the model. 



III. THE MODEL AND THE RESULTS 

Let us first of all illustrate the solvable model to which we applied the enlarged version 
of IRPA. 

It consists of three levels, 0, 1 and 2, with energies eo, £i and respectively. Let 2f2 be 
the degeneracy of each level and N = 2Q the total number of fermions in the system. 
We define the operators: 

^■y = E a lm a jm 1 (28) 

m 

where the indices i and j denote one of the three levels and the index m runs over the 2fl 
substates of each of them. The operators K satisfy the following commutation relations: 

[Kij, K kl ] = 5 jk Ku - 5 u K kj . (29) 

They are therefore the generators of the 27(3) algebra. The algebra becomes SU(3) if we 
consider the additional relation, 

N = J £K ii , (30) 

i 

that fixes the total number of particles. 

We introduce the hamiltonian for our system as follows: 

+V 2 E (K l0 K jk + K kj K 0l ) + ^3 E K » K ki ■ (31) 
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The terms with the Vo and V\ strengths describe, respectively, the phph and pphh parts of 
the interaction. The term with the V 2 strength is related to the ppph part, while the last 
term represents the pppp part. In standard RPA the only two-body terms of H that enter 
in the expressions for the matrices A and B of the equations of motion (^) are those with the 
strengths V and Vx, i.e. the ph two-body terms. In the I RPA approach [ |TT1 , P~2"| the ground 
state that is actually used in the calculations is correlated; so the single particle occupation 
numbers aren't strictly 1 for hole states and for particle states, as in \HF). In this case 
also the pppp term enters in the expressions for the matrices (actually only in the matrix 
A). In the I RPA approach with the enlarged configuration space all the terms contribute. 

The exact results for the system can be obtained either by using the SU(3) symmetry 
of the model or by diagonalizing the hamiltonian (|3l| ) in the complete set of states: 



\n x n 2 ) =C(K 10 r(K 20 ) n2 \0) , (32) 

where |0) denotes the state in which all the particles are in the level 0, n\ and n 2 are the 
numbers of particles in the levels 1 and 2, respectively, and C represents a normalization 
factor. 

With the same set of parameters chosen for the exact calculation, after having performed 
a standard RPA calculation, we solved the equations of motion both in the I RPA approach 



of pd]JT2| and in the new approach, with the enlarged configuration space. In the IRPA 
case the operators Ql are defined as linear combinations of ph (iO, with: i ^ 0) and hp (Oi, 
with i 7^ 0) configurations, as in eq. ([!]): 

Ql = J2(X»K i0 - Y»K 0i ) , (33) 

i 

while in the enlarged calculation they are defined as in eq. (pi]) : 

Ql = YX x ti k *i - Y "Ai) > ( 34 ) 

i>j 



where: 



k * = (^< /2i ^ • (35) 
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With the definition (|35|) we get, for the excited states of the system, the same orthonormality 
conditions as given in eq.(§). 

Note that in (|3lD the indices i and j run over all the 3 single particle levels of the model. 



In both cases ( p3|) and (|3l) we have solved the non linear problem of eq.s(^) by means of an 
iterative procedure. We fixed the number of particles N equal to 10. In this case the number 
of exact eigenstates of the hamiltonian is 66. The RPA and IRPA calculations will give two 
excited states, since their configuration space is composed only by the two configurations 
(1,0) and (2,0). The enlarged IRPA will give three states, since its configuration space is 
made by the three configurations (1,0), (2, 1) and (2,0). 

We tested various values for the four parameters V , Vi, V2 and V3 and for the energies 
of the levels, the results being qualitatively the same. In Fig. 1 we show one case, where: 

e = 0, e x = e, e 2 = 2.5e , (36) 

V = -X, V 1 = X , V2 = ^ y 3 = ^- (37) 

Both the e and x parameters have the units of an energy. In the figure the excitation energies 
are represented versus the increasing values of the strength \. Dashed lines refer to exact 
values of energies. Among all the 66 exact eigenvalues the two represented ones are those 
with energies equal to 1 and 2.5 at % = 0; i.e. those which correspond to the two RPA and 
IRPA excited states. Dotted lines refer to RPA and dot-dashed lines to IRPA values. The 
three values corresponding to the enlarged IRPA approach are represented by full lines. 

The collapse point of RPA, where its first excitation energy becomes imaginary, appears 
at x — 0.024. We observe that both the IRPA and the enlarged IRPA calculations push 
the collapse point towards greater values of the strength parameter x; so, in this regard, 
both methods improve the RPA results. Moreover it can be seen that the two exact values 
are better approximated by the first and the third states obtained in the enlarged IRPA 
approach, than by the two IRPA states. This is especially evident for the higher state. 

It is interesting to focus the attention on the presence of the additional state that the 
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enlarged IRPA gives, with respect to RPA and IRPA. Actually this state doesn't corre- 
spond to any of the found exact states. This fact seems to indicate that it is a spurious state. 
On the other hand its energy is not zero or very small, as it happens normally for spurious 
states. Its energy always starts, when \ starts from zero, from the energy difference between 
the levels 1 and 2, and so it depends on how we fix the values e\ and 62- This would mean 
that, if we applied our approach in a realistic calculation, the spurious states that would 
appear wouldn't be easily recognized and eliminated, not having in principle zero or very 
small energies. This could cause problems in the interpretation of the calculated spectrum 
of excitations. A similar situation is encountered in RPA at finite temperature and was 
also found in ||. Let us look at the transition probabilities related to this state. Figure 
2 shows the transition probabilities related to the obtained states, for RPA and EIRPA 
calculations, for different x strengths. Let's observe in the figure the transition probability 
related to the spurios state, 



where |^ ) is the ground state, ) is the spurious state and F the one body operator QI8D 
with all the / a /?'s equal to one. We Cclll S66 ttlclt Psp is very small, with respect to the other 
two transition probabilities, only when x is f ar from the collapse point; when x approaches 
the collapse point the transition probability ( |3~8D becomes appreciable (see the case x — 0.04 
in the figure). The same trend is found for other sets of parameters. 

This means that in the evaluation of any physical quantity, in a realistic calculation, 
the existence of spurious states would have some influence and it would be important to 
recognize and eliminate them from the calculation. The problem of how to recognize them 
is still open, as they don't have in general small energies and/or small transition probabilities. 

We present now, in table 1, the results obtained for the EWSR, in IRPA and in enlarged 
IRPA cases. These results refer again to the choice fl36|) and (|37| ) for the parameters. The 
violation of the EWSR, in the IRPA case, increases with the increasing values of the x 
parameter and is about of 30% for x — 0.03. In the enlarged IRPA case the EWSR is 




(38) 
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always exactly satisfied, as expected. 

This is an important achievement of the present method, because, as stressed in Section 
1, all the methods that have been proposed so far in order to go beyond the RPA, by 
avoiding the inconsistency of the Quasi Boson Approximation, always violate the EWSR 
identity. 

IV. CONCLUSIONS 

In conclusion, we have presented an extension of RPA which avoids the use of the 
Quasi Boson Approximation and, at variance with many other attempts made in the same 
direction, preserves exactly the EWSR. This is obtained as a generalization of a previously 
studied approach by enlarging the configuration space with respect to that commonly used, 
which contains only particle-hole elementary escitations. The approach has been tested on 
a 3-level solvable model. 

The authors gratefully thank P. Schuck and M. Sambataro for helpful discussions. 
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FIGURES 

FIG. 1. Excitation energies versus ^, with the (|36|)-(^7|) parameters. The energies in Y axis 
are expressed in units of e. 

FIG. 2. Transition probabilities for four values of j, with the (^)-(|37|) parameters. The 
energies in X axis are expressed in units of e. 
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TABLES 





l.h.s. IRPA 


l.h.s. Enl. IRPA 


r.h.s. 


0.012 


1.84957 


2.1877893606 


2.1877893605 


0.03 


0.89411 


1.3607871867 


1.3607871868 


0.04 


0.67858 


1.5861229231 


1.5861229230 



TABLE I. The l.h.s of EWSR in the IRPA and in the enlarged IRPA cases, together with 
the r.h.s., for different values of x/ e an d with the parameters (|33"|)-(|34|), 
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